# include "scene.h"
# include <stdio.h>
# include <malloc.h>
# include <math.h>

# ifdef _OPENMP_
# include <omp.h>
# endif

extern short conf_grid(struct grid_config *gc, struct grid *g);
extern short grid_gen(struct grid_config *gc, struct grid *g);

# ifdef SOLVE_ENERGY
extern short heater_param(struct grid *g);
# endif

extern void mem_alloc(FILE *lptr, struct arr_struct *arr);
extern void mem_distroy(FILE *lptr, struct arr_struct *arr);

extern short solve(struct grid *g, const struct arr_struct *depVar, const struct arr_struct *memFace, const struct arr_struct *memOut);

int main(int argc, char *argv[])
  {
   struct grid g;
   struct grid_config gc;
   struct arr_struct tnext, depVar, memFace[3], memOut;;
   unsigned short dim, tn, nCPU;

   printf("Program starts: in 'main()'...\n");

   g.lptr = stderr;
   if( (g.lptr = (LOG_FILE ? fopen(LOG_FILE, "w"):stderr)) == (FILE *)0) g.lptr = stderr;

#  ifdef TRACE_FILE
   g.trptr = stderr;
   if( (g.trptr = (TRACE_FILE ? fopen(TRACE_FILE, "w"):stderr)) == (FILE *)0 ) g.trptr = stderr;
#  endif

   g.nTRD = nCPU = 1;

#  ifdef _OPENMP_
   fprintf(g.lptr, "Parallel computation using OpenMP.\n");

/* Setting number of threads more than the number of available physical CPUs. */
   g.nTRD = (nCPU = omp_get_num_procs()) + 0;
   omp_set_num_threads(g.nTRD);

/* Enabling nested-threading. */
   omp_set_nested(1);
   if(!omp_get_nested())
       {
        fprintf(g.lptr, "ERROR:: Nested threads are not enabled, which is required.\nExiting...\n");
        return -1;
       }

   fprintf(g.lptr, "Number of thread(s) used for this computation: %d\n", g.nTRD);
#  endif

   fprintf(g.lptr, "Number of available CPU(s): %d\n", nCPU);

/* --- Grid generation--- */
   if(conf_grid(&gc, &g) < 0) return -1;

/* Allocating arrays to store the grid-POINTS. */
   g.z = (double *)malloc((g.nz + 2*NzG + 1)*sizeof(double));
   g.r = (double *)malloc((g.nr + 2*NrG + 1)*sizeof(double));
   g.th = (double *)malloc((g.nth + 2*NthG + 1)*sizeof(double));

   if(grid_gen(&gc, &g) < 0) return -1;

/* --- Grid generation done --- */

#  ifdef SOLVE_ENERGY
   fprintf(g.lptr, "NOTE: The energy equation will also be solved. Heater will be ON after %E seconds of the flow-start.\n", t_HEATER_ON);

/* Finding the heater's radius on the top surface. */
   if(heater_param(&g) < 0) return -1;
#  endif

/* Allocating memory. */

/* Allocating memory to store 3xp and 4x(u, v, w), and 1xT, while solving the Momentum and Energy equations. */
   depVar.n1 = g.nr + 2*NrG;
   depVar.n2 = g.nth + 2*NthG;
   depVar.n3 = g.nz + 2*NzG;
   depVar.nq = 15;
# ifdef SOLVE_ENERGY
   depVar.nq = 16;
# endif
   mem_alloc(g.lptr, &depVar);

   g.u = depVar.q[0];
   g.v = depVar.q[1];
   g.w = depVar.q[2];
   g.p = depVar.q[3];
#  ifdef SOLVE_ENERGY
   g.T = depVar.q[15];
#  endif

/* Allocating memory to store A, B and C's at the cell-faces */
   for(dim = 0; dim < 3; dim++)
           {
            memFace[dim].n1 = g.nr + (3-dim)/3;
            memFace[dim].n2 = g.nth + dim%2;
            memFace[dim].n3 = g.nz + (1+dim)/3;
            memFace[dim].nq = nRK;
            mem_alloc(g.lptr, &memFace[dim]);
           }

/* Allocating memory to store u,v,w,T at the outlet for tn_BACKUP number of time-steps.
   It is required for application of the Orlansky BC. */
   memOut.n1 = g.nr_gap;
   memOut.n2 = g.nth;
   memOut.n3 = zCells_BACKUP;
   memOut.nq = (DEPVAR_N-1)*tn_BACKUP;
   mem_alloc(g.lptr, &memOut);

   for(tn = 0; tn < tn_BACKUP; tn++)
     {
      g.outBC.u[tn] = memOut.q[tn];
      g.outBC.v[tn] = memOut.q[tn_BACKUP+tn];
      g.outBC.w[tn] = memOut.q[2*tn_BACKUP+tn];
#     ifdef SOLVE_ENERGY
      g.outBC.T[tn] = memOut.q[3*tn_BACKUP+tn];
#     endif
     }


/* Allocating memory to store dependent variables at outlet-cell. It will be used in application of the Orlansky BC. */
   tnext.n1 = DEPVAR_N-1; /* Corresponding to u, v, w, T */
   tnext.n2 = g.nr_gap;
   tnext.n3 = g.nth;
   tnext.nq = 1;
   mem_alloc(g.lptr, &tnext);
   g.outBC.tnext_ug = tnext.q[0][0];
   g.outBC.tnext_vg = tnext.q[0][1];
   g.outBC.tnext_wg = tnext.q[0][2];
#  ifdef SOLVE_ENERGY
   g.outBC.tnext_Tg = tnext.q[0][3];
#  endif

/* Memory allocation done. */

   solve(&g, &depVar, memFace, &memOut);

/* Freeing memory. */
   mem_distroy(g.lptr, &depVar);
   for(dim = 0; dim < 3; dim++) mem_distroy(g.lptr, &memFace[dim]);
   mem_distroy(g.lptr, &memOut);
   mem_distroy(g.lptr, &tnext);
//    mem_distroy(g.lptr, &coeff);

   free(g.r); free(g.th); free(g.z);


#  ifdef TRACE_FILE
   if(g.trptr != stderr) fclose(g.trptr);
#  endif

   fprintf(g.lptr, "...'main()' ends!\n");
   if(g.lptr != stderr) fclose(g.lptr);

   printf("...program ends.\n");

   return 0;
  }
